### Jonathan D. Klingler, Gary E. Hollibaugh, Jr., and Adam J. Ramey
### "Don't Know What You Got: A Bayesian Hierarchical Model of Neuroticism and Ideological Uncertainty
### Political Science Research and Methods
###
###
### Hierarchical Model Estimation Code
### Description: This file (one of eight) sets up the hierarchical model on a cluster node.
###
### Note 1: Run the Preprocessing_and_Summary.R file beforehand to generate the CCES2014.Rda file.
### Note 2: This file ends in SingleCoreXXX.R, and generates a file called DKWYGCoreXXX.Rda, where XXX is a number from 1-8.
### Note 3: Make sure the directories are set up properly and that all files are in this file's directory.
### Note 4: For ease of use, template script files called DKWYGCoreXXX.sh are included in the "Cluster Scripts" directory.
### Note 5: This file will start a process that will take several days to complete, so a cluster is recommended.


rm(list=ls())
ptm <- proc.time()

library(foreign)
library(rjags)
library(grDevices)
library(MASS)
library(random)



load("CCES2014.Rda")


# need to remove datapoints where questions were not asked
# and where personality not elicited
cces <- cces[!is.na(cces$self_emoti) & !is.na(cces$self_consc) & !is.na(cces$self_agree) 
             & !is.na(cces$self_extra) & !is.na(cces$self_openn) & !is.na(cces$CC421a) 
             & !is.na(cces$newsint) 
             & !(cces$CC334C == "Not Asked") 
             & !(cces$CC334D == "Not Asked") 
             & !(cces$CC334E == "Not Asked") 
             & !(cces$CC334F == "Not Asked") 
             & !(cces$CC334G == "Not Asked") 
             & !(cces$CC334K == "Not Asked") 
             & !(cces$CC334L == "Not Asked") 
             & !(cces$CC334M == "Not Asked") 
             & !(cces$CC334W == "Not Asked"),]

attach(cces)

##the dependent variable
obama <- CC334C
clinton <- CC334D
cruz <- CC334E
paul <- CC334F
bush <- CC334G
dem <- CC334K
rep <- CC334L
teaparty <- CC334M 
scotus <- CC334W

y <- cbind(obama,cruz,clinton,paul,bush,dem,rep,teaparty,scotus)
for (i in 1:ncol(y)) y[,i] <- as.numeric(y[,i])

y[y>=8] <- 0
y[is.na(y)] <- 0
y <- y+1

##covariates for decisiveness and the saliency intercepts
self_neuro <- 8 - self_emoti
self_openn <- (self_openn - 1)/6
self_consc <- (self_consc - 1)/6
self_extra <- (self_extra - 1)/6
self_agree <- (self_agree - 1)/6
self_neuro <- (self_neuro - 1)/6
polinterest <- (newsint=="Most of the time")*1
unkinterest <- (newsint=="Don't know")*1
faminc_recode <- faminc
faminc_recode[which(faminc_recode == "$150,000 - $199,999")] <- levels(faminc_recode)[17]
faminc_recode[which(faminc_recode == "$200,000 - $249,999")] <- levels(faminc_recode)[17]
faminc_recode[which(faminc_recode == "$250,000 - $349,999")] <- levels(faminc_recode)[17]
faminc_recode[which(faminc_recode == "$350,000 - $499,999")] <- levels(faminc_recode)[17]
faminc_recode[which(faminc_recode == "$500,000 or more")] <- levels(faminc_recode)[17]
faminc_recode[which(faminc_recode == "$250,000 or more ")] <- levels(faminc_recode)[17]
income <- factor(faminc_recode)
income_notsay <- (faminc_recode == "Prefer not to say")*1
age <- 2014 - birthyr
age2 <- (age^2)/100
black <- (race == "Black")*1
hisp <- (race == "Hispanic")*1
other <- (race == "Asian" | race == "Native American" | race == "Mixed" | race == "Other" | race == "Middle Eastern")*1
fulltime <- (employ == "Full-time")*1
parttime <- (employ == "Part-time")*1
unemploy <- (employ == "Unemployed")*1
retired <- (employ == "Retired")*1
female <- I(gender == "Female")*1
education <- as.numeric(educ)

X <- cbind(1,
           self_openn,
           self_consc,
           self_extra,
           self_agree,
           self_neuro,
           female,
           age,
           age2,
           black,
           hisp,
           other,
           education,
           polinterest,
           unkinterest,
           income,
           income_notsay,
           fulltime,
           parttime,
           unemploy,
           retired)

##covariate for saliency regression (equal to 1 if respondent and stimulus are the same party)
Xs <- cbind((CC421a=="Democrat")*1, (CC421a=="Republican")*1, 
            (CC421a=="Democrat")*1, (CC421a=="Republican")*1,
            (CC421a=="Republican")*1, (CC421a=="Democrat")*1,
            (CC421a=="Republican")*1, (CC421a=="Republican")*1,
            (CC421a=="Republican")*1)

forJags <- list(y=y, x=X, Xs=Xs)

forJags$J <- nrow(y)
forJags$K <- ncol(y)
forJags$R <- ncol(X)
forJags$nCat <- length(table(y))

## priors
forJags$b0 <- rep(0,ncol(X))
forJags$B0 <- diag(.25,ncol(X))

## initial values
inits <- function(){list(beta_e=rep(0,ncol(X)),
                   beta_d=rep(0,ncol(X)),
                   beta_a=rep(0,ncol(X)),
                   beta_s=0,
                   sig_e=1,
                   sig_a=1,
                   tau0=seq(1,2,length=(forJags$nCat-3)))}


n.tune <- 10000
n.burnin <- 40000
n.saved <- 50000
n.thin <- 1


set.seed(1)
temp.model <- jags.model(file="MLV2.bug",
                         inits=inits,
                         data=forJags,
                         n.chains = 1,
                         n.adapt = n.tune)
update(temp.model, n.iter = n.burnin)
dkwyg.1 <- coda.samples(temp.model,
                        variable.names=c("beta_s","tau","beta_d","beta_e",
                                         "beta_a","beta_b","gamma"),
                        n.iter = n.saved,
                        thin = n.thin)

save(dkwyg.1, file="DKWYGCore1.Rda")
detach(cces)
proc.time() - ptm

